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ABSTRACT 

This is the hrst in a series of papers on the effects of dust on the forma¬ 
tion, propagation, and structure of nonlinear MHD waves and MHD shocks in 
weakly-ionized plasmas. We model the plasma as a system of 9 interacting fluids, 
consisting of the neutral gas, ions, electrons, and 6 grain fluids comprised of very 
small grains or PAHs and classical grains in different charge states. We formulate 
the governing equations for perpendicular shocks under approximations appropri¬ 
ate for dense molecular clouds, protostellar cores, and protoplanetary disks. We 
describe a code that obtains numerical solutions using a hnite difference method, 
and establish its accuracy by comparing numerical and exact solutions for special 
cases. 

Subject headings: diffusion — dust, extinction — ISM: clouds — ISM: magnetic 
helds — MHD — plasmas — shock waves 


1. Introduction 

Shock waves in weakly-ionized plasmas usually exhibit a multifluid structure, in which 
the charged and neutral particles behave as separate, interpenetrating fluids. If the shock 
speed is less than the speed at which the charged fluid communicates compressive distur¬ 
bances, the charged particles will be accelerated ahead of the neutrals in a “magnetic pre¬ 
cursor,” and the neutrals accelerated, compressed and heated downstream by collisions with 
streaming charged particles (Mullan 1971; Draine 1980, henceforth D80). The flow is contin¬ 
uous (a G- or G* shock) if the length scale for acceleration of the neutrals is large compared to 


the scale for cooling; otherwise the flow undergoes a discontinuous transition (a J shock) at 
a viscous subshock (D80; Chernoff 1987; Roberge & Draine 1990). The dynamics, chemistry, 
and emission from multifluid shock waves can differ profoundly from those of their single¬ 
fluid counterparts. Consequently, they have been the subject of numerous investigations (see 
Draine & McKee 1993 for the most recent comprehensive review). 

The energy and momentum of the “charged fluid” are usually dominated by the mag¬ 
netic held, the charged particles acting mainly to provide energy and momentum transfer 
between the held and the neutral huid. At low densities the coupling is dominated by ion- 
neutral scattering and shock waves can be modeled as a three-huid system; the neutral 
gas, the ion/magnetic held huid, and the huid of electrons (which move with the ions but 
need not have the same temperature; see D80). However when the density of the preshock 
medium exceeds n ~ lO^crn”^,^ the energy and momentum transfer is dominated by colli¬ 
sions between neutral particles and charged dust grains (Elmegreen 1979; D80; Nakano & 
Umebayashi 1980; Ciolek & Mouschovias 1993). If the grains are well coupled to the mag¬ 
netic held, they can be incorporated into the 3-huid picture as an additional species of ion. 
This is a reasonable approximation for PAHs and very small grains with radii Ogg ~ 10“^ cm 
(henceforth referred to collectively as “small grains”) at densities up to ~ 10® (Neufeld & 
Hollenbach 1994). If the small grains are present in dense gas with a fractional abundance 
not much less than ~ 10“^ (the abundance inferred for the dihuse ISM, see Li & Draine 2001 
and references therein), they will dominate the held-neutral coupling and a 3-huid model of 
the dynamics is appropriate (e.g., Kaufman & Neufeld 1996a,b). 

However if small grains are depleted in dense gas, as seems likely, then scattering by 
large grains with Og > 0.1 pm will dominate the held-neutral coupling. Because the large 
grains are only partially tied to the magnetic held for n > 10^ cm“^, they must be treated as 
an additional huid (Chernoh 1985; Havnes, Hartquist & Pilipp 1987; Ciolek & Mouschovias 
1993). The problem is further complicated by the fact that large grains have diherent charge 
states, a continuous spectrum of sizes, and masses many orders of magnitude larger than 
those of the ions. Draine (D80) developed an approximate description of the grain dynamics 
by treating the grains as test particles and neglecting their inertia; in this approximation 
each grain drifts through the neutrals with a speed that can be determined by balancing the 
gas drag and electromagnetic forces. Chernoh (1985) treated the grains as a separate huid, 
included their inertia, and demonstrated that the grains can undergo a J shock. Pilipp, 
Hartquist & Havnes (1990) and Pilipp & Hartquist (1994) carried out numerical simulations 
of steady 4-huid shocks, emphasizing the importance of grain charging processes and the 


^We use n to denote the density of H nuclei in all forms. 



dynamical effects of the grain charge and current. Wardle (1998, see also Wardle & Ng 
1999) has formulated equations of motion valid for steady multifluid shocks with an arbitrary 
number of grain fluids, in a form that is valid for very large densities, and has discussed 
optimal procedures for representing the spectrum of large grain sizes by a single, effective 
size. 


This is the first in a series of papers on the effects of dust on the formation, propagation, 
and structure of nonlinear MHD waves and shocks in dense molecular clouds, bipolar out¬ 
flows, and protoplanetary disks. Here we describe a dynamical model that incorporates and 
extends previous work in several ways: (i) We allow for time-dependent flow. Prior work on 
time-dependent, multifluid shocks (Toth 1994; Smith & Mac Low 1997; Mac Low & Smith 
1997; Stone 1997; Chieze, Pineau des Forets & Flower 1998) has been based on the 3-fluid 
model. We are interested in time-dependent solutions because they permit us to study the 
formation of magnetic precursors, because of their possible relevance to bipolar outflows 
(Flower & Pineau des Forets 1999), and because time-dependent simulations are a robust 
way to hnd steady solutions (e.g., steady J shocks) that are otherwise hard to compute, (ii) 
We include the inertia of large grains, which has been neglected in all prior studies except 
Chernoff’s (1985). Neglecting the inertia is a good approximation when the gas drag time 
for large grains is small compared to all other dynamical time scales. While this is usually 
appropriate for the dynamical time of the neutrals, it is often a poor approximation for the 
charged fluid (see §2.3.3). A realistic treatment of grain dynamics is necessary for studies on 
the formation of magnetic precursors: Draine & McKee (1993) pointed out that magnetic 
precursors could not form (and C shocks therefore could not exist) for shocks faster than 
^ 20 km s“^ if large grains are able to “load” the magnetic field lines efficiently. The code 
we describe was designed in part to address this issue, (iii) We include small grains as a 
separate fluid. This permits us to model grain charging realistically and to examine vari¬ 
ous scenarios for small-grain depletion in dense gas. (iv) We include charge fluctuations by 
treating the grains in different charge states (3 each for small- and large grains) as separate 
fluids. Charge fluctuations exchange mass, momentum, and energy between the charged and 
neutral grain fluids, with implications for the loading of magnetic held lines. Nevertheless, 
the scope of our calculations is limited in some important respects: We do not attempt to 
describe the continuous size distribution of the large grains, which are assumed to have uni¬ 
form radii. We have also restricted our initial investigations to one-dimensional geometry so 
that, for example, time-dependent simulations of Wardle instabilities with dust are beyond 
our present capabilities. 

The plan of this paper is as follows. In § 2 we give the equations of motion and discuss 
the conditions where various approximations are likely to be valid. In § 3 we describe 
our numerical methods and summarize a suite of benchmark calculations that establish the 
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accuracy of our code. Our results are summarized in § 4. 


2. Formulation 


2.1. Modeling Assumptions 


We consider plane-parallel disturbances propagating along the +z direction of a Carte¬ 
sian frame, {x,y,z), defined so that the undisturbed plasma a.t z = -|-cx) is at rest. Fluid 
variables are assumed to depend only on and t. The plasma is threaded by a magnetic 
held of the form B = Bx{z,t)x. That is we restrict the discussion to perpendicular shocks; 
oblique shocks will be described elsewhere. Fluid motions in the x-direction (i.e., along the 
held) are uncoupled from those in the y and directions and are ignored. 

We model the plasma as a system of 9 interacting huids. Fluid variables are labeled 
by a subscript f, where f = n, i, e, sg“, sg°, sg"*", g~, g°, and g+ denote the neutral huid, 
ions, electrons, small grains with charge —e, 0, or +e, and large grains with charge —e, 0, 
or -|-e, respectively. The particles of huid f have mean mass mf, number density Uf, mass 
density pf, how velocity Vf, pressure Pf, and internal energy density Pf, where Pf includes 
the contributions of thermal motions, rotational and vibrational excitation, dissociation, etc. 


The neutral and ion huids are composed of W and W distinct chemical species, respec¬ 
tively, so that 




Un — ^ ^ ^na 
a=l 


( 1 ) 


and similarly for the ions. Realistic simulations require chemical networks elaborate enough 
to predict the abundances of important coolants (primarily H 2 O, CO, and H 2 in dense 
plasmas; see Neufeld & Kaufman 1993; Neufeld, Lepp & Melnick 1995) and to model changes 
in the fractional ionization (which can profoundly ahect the hydrodynamics; see Flower, 
Pineau des Forets & Hartquist 1985). In subsequent papers we include chemical networks 
tailored to the particular applications of interest to us. Here we are interested mainly in 
establishing the accuracy of our numerical code. The benchmark calculations we present 
were obtained using a simple analytic expression for the radiative cooling rate (see §3); 
the chemistry of coolants is therefore not required. To model the ionization we adopt a 
drastically simplihed “chemistry” with just 4 neutral species and 2 ions. The neutral fluid 
is composed of molecular hydrogen with uhj = 0.5n, atomic helium with riue = O.ln, plus 
a generic heavy atom, a, and molecule, m. The ion fluid is composed of the atomic ion a"*" 
and molecular ion m+. Since the dominant ions in the dense astrophysical plasmas relevant 
to star formation and protostellar disks are HCO^, Na^, and Mg’*", we take the ion mass to 
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be ma+ = min+ = rrii = 25amu. The chemical abundances are determined by a network of 
20 gas-phase and grain surface reactions (Appendix A). The fractional ionization predicted 
by our simple network is in good agreement with more elaborate models for n > 10 ^ cm“^ 
(Fig. 1). 

The large and small grains are assumed to have uniform radii Og and Ogg, where Og ~ 
10“^ cm and Ogg ~ 10“^ cm are adjustable parameters. Modifying our code to include a range 
of grain sizes would be straightforward but is beyond our current computational resources. 
The total abundances of large and small grains, summed over all charge states, are assumed 
constant. 


2.2. Governing Equations 

2.2.1. Neutral Particles 


The mass, momentum, and energy equations for the nentral flnid are 

dpn d 
dt + dz 
d d 

(Pn^nz) T (Pii^n T T Qvisc) Fnzi 


dt 


dz 

d d N rz 

yPiNny) T yPrNny'Vnz) ~ Fnyj 


and 


m 


rPn'i^n T Fn ) -|- 


dz 


rPn'i^n T Fn P Pn F Q\ 


— T’n'l’n + Gn — Aj; 


( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 


respectively, where = Vm'^ + Vny'^; the “sonrce terms” Sq, and Gn are the net rates 
per volume at which mass, momentum and thermal energy are added to the nentrals by 
interactions with the other flnids, and An is the rate of radiative cooling. The benchmark 
tests presented in §3 were carried ont for an idealized fluid with An prescribed analytically and 
En = Pn/ (7 — 1), with constant 7 . Fnture papers will incorporate a realistic calculation of 
radiative cooling by CO, H 2 O, and H 2 (Neufeld & Kanfman 1993; Neufeld, Lepp & Melnick 
1995) and the internal energy of H 2 acconnting for its time-dependent excitation and H 2 
dissociation. Our model includes artihcial viscosity to smear out jump fronts in the nentral 
flow. We use the von Nenmann-Richtmyer prescription for the psendo-viscous stress. 


q{Az)‘^ Pn \dvn/dz\^ if dvn/dz < 0 , 


Qy'\sc 


( 6 ) 


0 


otherwise 
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(see Richtmyer & Morton 1967), where I\z is the grid spacing and q is an adjustable dimen¬ 
sionless parameter. If present, a jump front will be smeared over ~ q grid cells. 


2.2.2. Electrodynamics 

We assume that the magnetic field is frozen into the fluid of ions and electrons. (The 
validity of this approximation is discussed in §2.3.) In the frozen held approximation, the 
electrons and ions move with a common velocity, 

W = (7) 

the electric held satishes 

E^-—xB, (8) 

c 

and the Maxwell equations reduce to the induction equation, 

plus the requirement of macroscopic charge neutrality, 

ne = ni + ng+ 4- nsg+ — n^- — Ugg-. (10) 


2.2.3. Electrons 

We do not solve mass or momentum equations for the electrons because Uf. and are 
functions of the other huid variables in the frozen-held approximation (cf. eqs. [7] and [10]). 
The electron energy equation is omitted in favor of Chernoh’s (1987) approximation, 

Ge = Ae, (11) 

which says that the electron huid rapidly reaches thermal balance due to its small heat ca¬ 
pacity. The shock structure depends on the electron temperature via the rate coefficients 
for electron impact excitation of molecules (which contribute to the radiative cooling) and 
the rate coefficients for recombination and grain charging (which inhuence the ionization; 
see Appendix A). The solutions presented in §3 are insensitive to Te because the cooling is 
prescribed analytically and the how time is much less than the timescale for changes in the 
ionization. We therefore adopt the (unrealistic) approximation Tg = O.lSTj for expediency; 
this approximation is based on the results of Draine, Roberge, & Dalgarno (1983, here¬ 
after DRD83). Future papers on realistic simulations will incorporate a detailed description 
of electron heating and cooling including the atomic and molecular processes described in 
DRD83. 
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2.2.4- Ions 


We solve the mass equation for the ions, 

^ ^ (PiVi.) = -^n, ( 12 ) 

but not the momentum or energy equations. We calculate the momentum of the ion fluid by 
(i) neglecting the inertia of the ions and electrons and (ii) requiring the total electric current 
to be consistent with the magnetic held gradient. We assume that the ion and electron 
velocities are determined solely by the balance between electromagnetic forces and drag due 
to elastic collisions with the neutrals. Then requirement (i) implies 




0 = -eUe [E -\ - xB] + —5- (Un - u, 


Pe 




(13) 


and 


0 = enJE + —xB) + {v^ - v^) , 


Pi 


(14) 


where rd is the drag time for species i due to elastic scattering with species j (a function of 
the huid variables, see DRD83).^ Adding expressions (13) and (14) gives 


6 Pi 

e (Ui - He) E + - {UiVi - UeVe) XB + -^ {v^- Ui) 
C 


i,n 


0, 


(15) 


where we have neglected the drag force on the electron huid compared to the drag force on 
the ions.^ The quantities Ui — and n\V\ — UeVg can be written in terms of the densities 
and velocities of the other charged particles using the neutrality condition, (10), and the 
dehnition of the total current. 


j = {UiVi — UeVe + ng+Vg+ — Ug-Vg- + nsg+Vsg+ — Usg-Vsg-) e. 

Making these substitutions in eq. (15) we hnd, after a little algebra, that 
Pi {'VxB)xB en. 


(16) 




Un - VP = 


All 


eug- 


[ivg,-vOxB]-^[i 


PT) — 

[(Usg- - Ui) xB] H-^ [(Usg- - Ui) xB ], 


u 


sg^ 


- -Wi) xB] 


(17) 


^Notice that equations (13) and (14) would yield the nonsensical solution Vg = Vi = if one naively 
inserted expression (8) for the electric held. The meaning of this result is that the drag force must vanish 
in the frozen-held approximation, where gas drag is taken to be negligible compared with electromagnetic 
forces. 

^The ratio of the drag forces is of order {me/mn)crenVe /^ 10”^, where (Jen is the geometrical cross 
section, Vg is the electron thermal velocity, and (cryj^^ is the ion-neutral momentum transfer rate coefficient 
(DRD83). 
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where we used equation (8) to express the electric held in terms of Vi and Ampere’s Law 
to write the current in terms of VX-B. Equation (17) says that the electromagnetic and 
drag forces on the combined ion+electron huid are in balance, as they must be when ion and 
electron inertia are neglected. Equation (17) is solved to hnd V[{z,t) in terms of the other 
huid variables. 

The ion energy equation is omitted in favor of the approximation 

/CbT; ^ /CeTn + I^»i - (18) 

(Chernoh 1987), which says that the ion huid is in thermal balance with heating and cooling 
due to ion-neutral scattering. 


2.2.5. Small Grains 


Our code solves the equations for mass conservation. 


dpsg- , / \ Q 

^ dz ^ ’ 


dt 


(19) 


dpsgo , d 


dt dz 5'sgO, 


and 


dp, 


{Psg+Vsg+z) — *S'sg+. 


sg^ 


( 20 ) 


( 21 ) 


The source terms on the r.h.s. of equations (19)-(21) describe the transfer of grains from 
one charge state to another due to the capture of ions and electrons, and are computed by 
integrating the appropriate cross sections (Draine & Sutin 1987) over the distribution of 
relative velocities for charged particles and drifting grains (see Appendix A). We neglect the 
inertia of the small grains and compute their velocities algebraically by requiring the total 
force on each small-grain huid to vanish. Thus we require 


0 = Psg- {v 


Sg 


"t^i) xr 2 sg- -|- 


r 


'^sg-y 


0 

cf V- r t^sgU 


sg ,n 


r 


sg ,1 


V,e 


( 22 ) 


0 = 


PsgO 
'sg®,n 


(^n - ^sgo)-^^sgO 


^cf 

^gO,i 


PsgO 

T-cf 

hgO,e 


UsgO -|- 


Psg- 

sg ,1 


n. 


sg- 


+ 


Psg+ 


T, 


cf 


V 


sg-f 


sg-c,e 


(23) 


Ui) X L2sg+ T 


Psg+ 


(Un - nsg+) 


sg“^,n 


PsgO 

cf "^SgO 

AgO.i 


Psg+ 


r 


cf 


V 


SgT, 


sg^,e 


(24) 


and 


0 — Psg+ ('*^sg 



9 


where 17 f (= QB/rrifc) is the cyclotron frequency for particles of type f with charge Q 
and Tj^J is the mean time for a small grain to change its affiliation from fluid i to fluid 
j by capturing an ion or electron. As noted in previous studies (e.g., Pilipp et ah 1987; 
Ciolek & Mouschovias 1993), the capture of charged particles by the grains transfers mass 
and momentum between the different grain subfluids. Our code includes these effects via 
the charge fluctuation terms in the momentum equations and the source terms in the mass 
equations. The velocity dispersion of the small grains does not affect other properties of the 
flow; consequently, we do not solve the small-grain energy equation. 


2.2.6. Large Grains 

We neglect the velocity dispersion of the large grains and omit the energy equation. We 
include the large-grain inertia and solve the equations for mass conservation. 


dt + 5;, 
dt dz 

(9pg+ d , ^ r, 


(25) 

(26) 
(27) 


and momentum conservation. 


m 


d 


(Pg W y) A (Pg y'^g z) Pg (w 2: 'Viz) ^g + d ( 


dz 


r _ 
g ,11 


Vny Vg-y) 


P^V _ + P^Ly 0 
cf ^g y ^ rf 


T, 


V,e 


(28) 


d 


d 


(PgOVgOy) + (PgOVgOyVgOz) . {Vay VgOy) ^ , V g- 


Pg- 


Pg° 


yrg'^'^g'^y/ I Krg'^'^g'^y'^g 


V,n 


r 


cf y 


cf 

V,e 


.^cf "^g^y 
Ao ,■ 


/5g+ 

V„Oy -\ - 


T, 


cf ^g^y^ 


g~^,e 


(29) 


d d 

^ (Pg+'*^g+?;) + ^ (Pg"W+ 2 /W+^) 


~ Pg+ {Vg+z Viz) f^g + ^ (Wy Vg+y) 


g^,n 


PgO 

H— ^v^o 


Pg+ 


_cf '^g'^y _cf '^g^y 

V,i V,e 


( 30 ) 
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m 


d 


{Pg-'^g-z) + Q^{pg~'^g-z) ~ Pg- {'^g-y 


r _ 
g ,n 


0 

cf Z “ cf g 


r 


■ ,1 


v.» 


(31) 


and 


<9 . X <9 / 2 ^ 

(Pgongo,) + — (pgongoj 






Pg^ f \ Pg~ PgP 

“d ~ '^g°z) + cf '^g~ ^ ~ 'Z^'^g°^ 

g°,n ^g-,i ^gO,e 

Pg° Pg+ 

-VsOz H—■ 


cf '^g'^z I cf '^g’^2’ 

V,i V,e 


<9 / . d , 2 \ 

(Pg+ng+z) + — (pg+Wg+j,) 


Pg+ ("^g+j/ ^g 9“ d ^g+z) 


T 


I '8 o, 

“T „r VnOr 


_cf '^g" 
VO 1 


Pg+ 

-^t^g+z, 

V,e 


where = Ifip^-I = I fin 


2.2.1. Chemistry 


(32) 


(33) 


The ion and nentral mass density depend on the chemistry, which is calculated in parallel 
with the hydrodynamics. Our code solves a set of chemical equations for the neutrals, 


and ions. 


Qji Q 

A- ^ (n y ) = c —V 

r\, ~ r^ y PrVOi^YVZ) 

at oz 

OtI' 3 

-L _ (n- V- ) = C- -V- 

\ ^ y^\Oi^\Z) ^lOL ^lOLI 


(34) 


(35) 


dt dz 

where Cfa and Vfa are the net rates per volume at which species a in fluid f is created and 
destroyed, respectively. 


2.3. Validity and Relevance 

2.S.L Frozen Field Approximation 


The frozen held approximation is valid when Tj 3> 1, where 

eB il_3 

Ti =-^36 —^ 

ns 


m\C 


( 36 ) 
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is the ion Hall parameter, B = i ?_3 mG is the magnetic held, and n = lO^nscm"^. The 
line-of-sight magnetic held components inferred from Zeeman splitting in thermally-excited 
transitions (summarized by Crutcher 1999) imply that Tj 3> 1 for n ~ lO^cm”^ but observa¬ 
tions of denser regions are rare. Extrapolating the thermal Zeeman data to higher densities 
using the observed correlation B oc (with k, ~ 0.5, see Crutcher 1999; Basu 2000) would 
imply that the frozen-held approximation is valid for n > 10®cm“^. The extrapolated held 
values are consistent with the helds estimated in HoO masers with n > 10®cm“^ but the 

^ rs-/ 

held strengths inferred from maser observations are controversial (Sarma, Troland & Romney 
2001 and references therein). 


2.3.2. Inertia of Ions and Small Grains 


Neglecting the inertia of ions and small grains is an excellent approximation for steady 
shocks. Let Tdrag(sg) be the length scale for gas drag to decelerate the small grains, Tneut the 
scale for compression of the neutrals, and Tfieid the scale for compression of the magnetic held. 
Figure 2 compares these and other length scales for a shock with speed Vs = 10 km s“^ and 
various preshock densities. The dynamical scales were estimated using dimensional analysis 
of the equations of motion for a steady shock (DRD83; see eq. [3]-[4]): 


L 


neut 


PnO Vs 
Fn 


(37) 


and 


Lfield — 


SttFi ’ 


(38) 


where the subscript “0” denotes preshock quantities. We took the preshock magnetic held 
to be 


Bq = b (no/cm pG 


(39) 


with b = 0.3. We included the drag due to electron-neutral and ion-neutral scattering 
(DRD83, eq. [11]-[14]) as well as the gas drag on small and large grains (D80, eq. [40]). 
We computed Lfiew at a typical point in the magnetic precursor where the neutral particles 
are unperturbed and the magnetic held has been compressed by a factor of 2 (see D80, 
eq. [19]). For Lneut we took a representative downstream point, with unperturbed neutrals 
and the magnetic held compressed by a factor of M^a, the neutral Alfven Mach number (see 
Kaufman & Neufeld 1996a and references therein). Preshock chemical and grain abundances 
were taken from Figure 1 and the drift velocities of large and small grains were computed 
by neglecting their inertia for expediency. We set 


-hdrag(sg) '^sg-z T'sg-,n'> 


(40) 



12 


with the grain velocity evaluated in a frame where the shock is steady and the drag time 
computed from the velocity-dependent drag force on a drifting sphere (Draine & Salpeter 
1979). /^From Figure 2 and analogous calculations for other shock speeds we conclude that 
the inertia of small grains (and hence the ions) can be neglected in steady, nondissociative 
shocks with Vs < 30 km s“^. 

In time-dependent shocks, the inertia of light particles is a subtler issue. The formation 
of a magnetic precursor depends on the speed at which compressive disturbances are com¬ 
municated to the charged fluid, and this “signal speed” is modihed when the inertia of light 
particles is neglected. As a simple illustration, consider the propagation of small-amplitude 
waves in a plasma of ions, electrons, and neutrals. In Figure 3 we plot the phase velocity of 
waves propagating perpendicular to B, as predicted by the exact dispersion relation for a 
cold plasma (solid curves) and the dispersion relation with ion inertia neglected (open cir¬ 
cles). The exact dispersion relation predicts the well-known result (Kulsrud & Pearce 1969) 
that compressive waves propagate in the ion-electron fluid for wavelengths 

A < Ai = dTTOims'h'^n (41) 

and in the bulk plasma for 

A > A 2 = Ai a/ Pn/ 4pi, (42) 

where U ims = B/ is the ion magnetosound speed for a cold plasma. When the ion 

inertia is neglected, the short-wavelength mode disappears and the fastest “signal speed” for 
propagating waves drops from Uims to the bulk magnetosound speed, Ums ~ \/p-JPn- This 

seems to suggest that the “inertialess equations” would fail to admit C type solutions for 
Vs > Urns- However Smith and Mac Low (1997) have shown, by direct numerical integration 
of the inertialess equations, that a J shock with Vg > Vms can evolve into a C shock. We 
describe a similar result in §3. The development of a magnetic precursor in both examples is 
initiated by rapid ion motions with ~ Vjms- Evidently the formation of a C shock does not 
require propagating waves with phase velocities larger than the shock speed. We speculate 
that the ions are accelerated by evanescent waves in the “forbidden zone” Ai < A < A 2 , 
which have long damping times (Fig. 4). If correct, our conjecture may have interesting 
implications for maximum speed of a C shock. This and related issues will be discussed 
elsewhere (Roberge & Ciolek, in preparation). 


2.3.3. Relevance of Large-Grain Inertia and Charge Fluctuations 

The drag length for large grains [Ldrag(g)] is shown in Figure 2. For n < 10®cm“^, the 
large grains are partially coupled to B and may affect the formation of C shocks by “loading” 
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the field lines (Draine & McKee 1993; see §1). Over this density range the drag length is 
comparable to -fifieid (Fig. 2b), so including the large-grain inertia is necessary. Length scales 
for the capture of charged particles by small (sg°-|-e, sg“-|-i) and large (g°-l-e, g“-l-i) grains 
are also plotted in Figure 2. The charge fluctuation times were computed using abundances 
from Figure 1 and rate coefficients from Appendix A. Because the shock structure is sensitive 
to changes in the ionization (Flower, Pineau des Forets & Hartquist 1985; Pilipp et ah 1990), 
it is necessary to include the charge fluctuations of small grains for n < 10®cm“^ (Fig. 2a). 
Figure 2b shows that large grains undergo many charge fluctuations as the magnetic field is 
compressed, partially coupling the neutral grains to the field lines. Consequently the charge 
fluctuations of large grains are also important. 


3. Solutions 

3.1. Initial and Boundary Conditions 

We model two types of problem with our code. In the first type (“free propagation 
problems”), the initial conditions are some disturbance, linear or nonlinear, specified on a 
finite computational domain, Q < z < Z. For free propagation problems we impose “non- 
reflecting boundary conditions” (discussed below) at z = 0 and z = Z. These boundary 
conditions allow the disturbance to propagate off the domain as if the boundaries were 
at infinity. In the other type of problem (“inflow problems”), the initial conditions are a 
uniform, static “undisturbed state” for 0 < z < Z and a steady flow is incident at the 
boundary z = 0. For inflow problems we specify the flow variables at z = 0 and impose 
nonrefiecting boundary conditions at z = Z. 

For inflow problems, the undisturbed state is specified by choosing the density Uq and 
magnetic field strength, Bq. In the exemplary calculations presented in this paper, the 
temperature of the undisturbed gas is not of interest. We assign the initial neutral, ion, and 
electron temperatures to reasonable but somewhat arbitrary values instead of, say, requiring 
the fluids to be in thermal balance. Once the gas density and temperature are specified, the 
initial densities nf^o of the various chemical species are determined by solving the chemical 
rate equations for a static system (see §2.2.7). To determine the abundances uniquely we 
must specify the fractional abundance of heavy nuclei {xg), and the initial abundances of large 
(xg^o), and small (xsg,o) grains. We typically set Xg ~ 10“®. The initial grain abundances are 
specified in practice by choosing the grain radii, Og and Ogg, the densities pg and psg of the 
grain materials, and the initial dust-to-gas mass ratios, Xg_o and Xsg.o- For diffuse regions 
of molecular clouds or protostellar cores, we typically set Og = 0.1 pm pg = 3.3gcm“^, and 
Xg,o = 10“^, yielding Xg = 2 x 10“^^. For the small grains we usually take Ogg = dA, 
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Psg = 2.2gcm“^, and Xsg,o = 5 x 10“^, yielding an abundance, a;sg,o = 2 x 10“'^, consistent 
with the constraints on small grains in diffuse gas (Puget & Leger 1989; Weingartner & 
Draine 2001; Li & Draine 2001). For models of shocks in dense gas, we adjust the grain 
radii and dust-to-gas ratios to account approximately for the possible depletion of small 
grains, coagulation of large grains, (e.g., Ossenkopf 1993), and depletion of large grains due 
to ambipolar diffusion during core formation (Ciolek & Mouschovias 1994, 1996). 


The actual problem we solve is a dimensionless one. The governing equations (§ 2) are 
put into nondimensional form by using as a unit of distance some specihed length L, with 
Z = TIL, TZ a constant. Velocities are normalized to Cno, the isothermal sound speed in 
the undisturbed neutral gas. (For free propagation problems, “undisturbed gas” means the 
gas outside the initial disturbance.) The unit of time is therefore L/Cno- Temperatures and 
densities are normalized to the temperature, and mass density, pnO; of the undisturbed 
neutral gas, and the unit of magnetic held strength is taken to he Bq. A consequence of our 
selection of units is that the dimensionless Alfven speed of the undisturbed gas, nA,no/C'no 
[= 5o/(47rpno)'/'C'no], appears as a characteristic parameter in the problem. 


The use of non-rehecting boundary conditions (also referred to as “transmissive” or “ra¬ 
diation” boundary conditions) allows the how to pass unimpeded through the computational 
boundaries, and prevents the rehection of spurious waves from the boundaries back into the 
computational domain. A comprehensive review of non-rehecting boundary conditions in 
numerical hydrodynamics is presented by Givoli (1991). We use the following formulation, 
due to Thompson (1987) and Vanajakshi, Thompson, & Black (1989). Examination of the 
governing equations (§2) reveals that they are all conservative transport equations of the 
form 


dU 


dz 


C = 0, 


(43) 


where 


Pn 


U = 


Pn'^ny 

Pn'^nz 


f PnUn + E, 


(44) 


is the (column) vector of conserved variables, J- is the corresponding vector hux of conserved 
quantities, and —C is the vector of source terms. The system (43) is transformed into the 
quasilinear form 


du' 


+ A- 


du^ 

dz 


+ C' 


0 , 


dt 


(45) 
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where 


W 


Pn 

Vny 

'l^nz 

En/Pn 


(46) 


is the vector of primitive variables and the Jacobian matrix, is readily determined. 

The eigenvalues of A are the characteristic speeds of the flow. Non-reflecting boundary 
conditions are implemented by imposing the set of equations 

+ Ek + Ik'C = 0, (47) 

at the appropriate boundary, where Ik is the left (row) eigenvector associated with the 
eigenvalue ^k and 

{ ^kh'dU-'/dz for outgoing waves, 

(48) 

0 for incoming waves. 

As demonstrated by Thompson (1987) and Vanajakshi et al. (1989), the conditions (48) for 
outgoing waves < 0 at z = 0, > 0 at z = Z) and incoming waves (.^fc > 0 at 2 ; = 0, 

^k < 0 at z = Z) eliminate spurious reflections and produce accurate numerical solutions at 
the artihcial boundaries. 


3.2. Numerical Method 

The dimensionless governing equations with nonreflective and/or inflow boundary con¬ 
ditions comprise a set of coupled, nonlinear, partial differential equations (PDEs). We solve 
the PDEs numerically using a technique known as the “method of lines” (e.g., Liskovets 
1965; Schiesser 1991). The computational domain is divided into N cells of uniform width 
= Z/N with cell boundaries at {zj = jAz, j = 0,1,..., iV}. Mean values of the variables 
are assigned to the center of each cell, i.e., to the points 

_ ^ {Zj ^ ^ 1 ^^ ^ _ 1 ^^^^ 

for j > 1. The mean value of variable / in cell j will be denoted by fj = f{zj). In the method 
of lines, spatial derivatives are approximated using hnite differences, thereby converting the 
PDEs into a system of ordinary differential equations (ODEs) of the form 

dt 




(50) 
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where dy/dt, y, and f are vectors of length WN, and W is the number of variables dehned 
at each cell. The ODEs are solved using standard stiff ODE solving subroutines readily 
available in the public domain. 


For cells inside the computational domain, spatial derivatives are approximated with 
the fourth-order differencing scheme 


dz 


^ -(l/12)(/j« - /,_2) + (2/3)(/,+i - 


A 2 


For cells adjacent to the boundaries, we use the three-point scheme 


dz 

dz 


(1/3)/2 + /i — (4/3)/j 


bO 


N 


Az 

(4/3)/bjv — fN — (l/3)/Af-i 
Az 


(51) 

(52) 

(53) 


where fbo and fbN refer to the values of variables at the boundaries z = 0 and zn = Z, 
respectively. Thus, derivatives at the boundaries depend only on variables inside the com¬ 
putational domain. They are accurate to 0{Az)'^, while the interior derivatives (eq. [51]) 
are accurate to 0{Az‘^). 

Advection equations, of the form df/dt = —d{fv)/dz, are discretized using the mono¬ 
tonic, upwind algorithm of van Leer (1979). Following this algorithm, advection terms are 
approximated as 


djfv) 

dz 


-I 3 


Az 


(54) 


where the flux at the right cell face is approximated by setting = {vj + Vj+i)/2, and 

Ij + if > 0 , 


fn = 


(55) 


fj+i - if < 0, 

and the flux at the left cell face is approximated by setting = (Tj_i -f- Vj)/2 and 

Jj-i + |(A/)j_i if > 0, 


_ _ 

Here we have introduced the “monotonic derivative,’ 


(56) 


if < 0. 


(A/), = 


Ai-Z.-i 

0 


if (/,+i-A(/,-Ai)>o, 

otherwise. 


(57) 
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For the boundary cell adjacent to = 0 we set = Vbo- If > 0, /^ = /m, otherwise 
fi is given by equation (56). Similarly, for the cell adjoining z = Z, vr = v^n- If vr > 0, 
equation (55) is used for if vr < 0, /r = fbN instead. 


3.3. Benchmark Tests 

To establish the accuracy of our numerical code, in this section we present the results 
of various test runs and compare them with known solutions. In all of the tests presented 
here, we set i7n = Pn/{l — 1 ), where 7 is the adiabatic index of the neutral gas, assumed 
constant. In models that include cooling, we use the analytical cooling function An dehned 
by Roberge & Draine (1990, eq. [3.18]). 


3.3.1. Shock Tube 

A standard test for numerical codes is the well-known shock tube problem (e.g.. Sod 
1978). In the absence of magnetic helds, the solution for this type of Riemann problem 
can be found exactly using analytic methods (Toro 1997 and references therein). The gas 
is adiabatic in the shock tnbe test, with 7 = 5/3. The magnetic held strength is made 
vanishingly small; in this limit, the ions, electrons, and grains move with the nentrals (dne 
to collisional forces) and do not appreciably affect the how, because the density of charged 
particles is mnch smaller than the nentral gas density. 

The shock tube test falls into the category of “free propagation” problems with nonre- 
hecting bonndary conditions at both edges of the spatial domain (see §3.1). The nentral gas 
is initially at rest. The initial density and temperature prohles are constant on both halves 
of the domain with a discontinnity at z = Zj^. In our test we chose = 12pno and 

T^{z) = 1.875Tno for 2 : < Z/2 and Pn{,z) = pno and Tn( 2 :) = T^q for 2 : > Z/2, where pno and 
Tno are scaling valnes that merely dehne the dimensionless nnits (see §3.1). 

For this test we turned oh the sonrce terms in the continnity eqnations for all hnid 
species. Because the gas how is taken to be adiabatic, the cooling fnnction for the nentral 
gas is set to zero. 

The results of the shock tube test for a dimensionless time ti > 0 are displayed in 
Figures Aa-Ad (circles). This particular model used 400 mesh points; for visnal clarity in each 
plot we display only every eighth data point. Figure 4a shows Pn{z,ti)/p^io as a function 
of z/Z. The ratio Vnz{z,ti)/Cno is exhibited in Figure Ab. The pressure Pn{z,ti) relative 
to Pn,o, and Ta{z,ti)/T^fi are shown in Figures 4c and Ad, respectively. For comparison. 
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also displayed in these figures are the exact Riemann solutions (solid curves). Examination 
reveals that our multifiuid code effectively reproduces the full structure of the shock tube 
problem, including the shock wave, rarefaction wave, and contact discontinuity which are 
present at this time. The r.m.s. relative errors of the numerically-computed variables are 
typically ;< 10 “^ — except at the oscillations near the base of the rarefaction wave and at 
the contact discontinuity. 

Figures 5a-d present numerical results for the same model run (circles) at a later time 
t 2 , which are again compared to the exact solution (solid curves). By this time, both the 
forward-moving shock wave and contact discontinuity, and also the backward propagating 
rarefaction wave, all present in Figures Aa-d, have moved “off stage”, that is, past the 
boundaries at z = Z and 0. We note that passage of the forward waves through z = Z 
does not produce any substantial numerical reflections at the boundary located there. Nor 
are there any spurious reflections generated by the backward rarefaction wave incident upon 
the opposite boundary at z = 0. This illustrates the efficacy of the non-reflecting boundary 
conditions (§ 3.1). 


3.3.2. Nonmagnetic Shock with Inflow Conditions 

The next test we present is an inflow problem. We again consider nonmagnetic, adiabatic 
(7 = 5/3) flow. Initially the neutral gas within the entire region [0, Z] is at rest with uniform 
density pno and temperature Tno. At t = 0 a high-velocity shock with v^z = SOCno is incident 
on the boundary at z = 0; this corresponds to a Mach number M = = 51.7, 

where Vg is the shock speed. The inflowing shocked material has density pn = 3.99pno, and 
temperature Tn = 835Tno. 

Figures 6a-d depict the profiles of the density, velocity, pressure, and temperature of 
the neutrals, normalized to the uniform initial state, at a time t > 0 . As before, the 
circles represent the solution from our multifiuid numerical code. The solid line is the exact 
Riemann solution. Comparison indicates that there is good agreement between the numerical 
and exact solutions. 


3.3.3. Counterpropagating Rarefaction Waves in a Magnetic, Flux-Frozen Fluid 

This test is used to verify that our numerical code has the correct MHD behavior in the 
limit of flux-freezing, that is, when the magnetic field is frozen into the bulk fluid. Effective 
freezing of the magnetic flux into the neutrals can occur when the characteristic minimum 
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flow time Tflow.min (= Az/vnz) is 3> T^tot) ^^6 drag time for neutrals colliding with all charged 
particles attached to the magnetic held. In this situation, information about the magnetic 
held is quickly communicated to the neutrals, through collisions with the plasma, on a time 
scale much smaller than the natural time scale associated with the how. When this happens, 
the magnetic held, plasma, and the neutrals can rapidly respond to one another, and thus 
move together as a single combined huid. In our code, the condition Tfiow,min/'rn tot 1 is 
readily satished by using a state with a large enough abundance of charged particles (since 
Tnf cx e.g., see eqs. [2c] and [10a] of Ciolek & Mouschovias 1993), and/or choosing 
a sufficiently large value Z for the size of the physical region being modeled (recall that 
A;^ = ZjN). 

The model we choose for this test is the counterpropagating MHD rarefaction how de¬ 
scribed in § 4 of Ryu & Jones (1995). This test consists of two oppositely-directed, adiabatic 
rarefaction waves in a magnetized huid. The initial state has a nonzero, discontinuous ve¬ 
locity prohle with Viiz{,z < ZjZ) = —Cno, and > ZjZ) = Cno- However, all of the 

other physical quantities are uniform in the domain [0, Z]: Pn{,z) = pno, ^n(^) = Jno, and 
Bx{z) = Bq. The Alfven speed in the undisturbed neutrals, nA,no, is equal to Cno for this 
model. Creation/destruction of charged species was not considered in the analysis of Ryu 
& Jones. Hence, to be able to make an accurate comparison to with their published model 
results, in this test we switched oh the source and sink terms that appear in the continuity 
equation for each charged species. 

Figures la-d exhibit the state of this test model at a time t > 0, when the forward- and 
backward-traveling rarefaction waves have opened a cavity in the density and magnetic held 
where mass and magnetic hux have howed away from the center. The curves with circles 
show (a) the density, {h) velocity, (c) pressure, and (d) the magnetic held prohles generated 
by our numerical code at this particular time. For comparison, also shown in each plot is 
the solution (solid lines) of Ryu & Jones (1995, their Fig. 3b), which they obtained both 
numerically and by using a specially-developed MHD Riemann solver. These hgures show 
that our multihuid code reproduces the limiting ideal MHD how. That the how in this 
model is indeed ehectively hux-frozen can be seen by comparing Figures 7a and 7d, which 
show that the mass and magnetic hux have behaved in the same way, preserving the local 
mass-to-hux ratio {pn!B^)/{pno/B q) = 1 throughout the evolution. There is, however, a 
numerical overshoot at the half-way point, = Z/2, which is the location of the contact 
discontinuity within the initial state. A similar overshoot was also found to occur in the 
numerical simulation of Ryu & Jones. 

Finally, we note that there are no spurious oscillations at either boundary in Figures 
7a-d. This conhrms that the non-reflecting condition imposed at the artihcial boundaries in 



our code is also properly dealing with the passage of MHD flows. 


3.3.4- Chernoff’s Solution for Adiabatic C Shocks 

We also test the effectiveness of our code in dealing with realistic, non-ideal MHD 
situations, such as when the magnetic field is not necessarily frozen into the neutrals, due 
to relatively small abundances of charged particles. To examine this regime, we modeled 
the formation of a C shock in a weakly-ionized plasma. As noted earlier, this problem for 
three-component plasmas (neutrals, ions, and electrons) has previously been considered in 
great detail for both steady (e.g., D80; DRD83) and also nonsteady (e.g., Toth 1994; Smith 
& Mac Low 1997) flows. 

For the case of time-independent flow, a detailed phase plane analysis of ‘energy- 
conserving’ (i.e., adiabatic) MHD shocks was presented by Chernoff (1987) for three-component 
plasmas. We use his results as the benchmark for this particular test. 

To compare our results to those of Chernoff, we used a model containing relative abun¬ 
dances of dust grains that are orders of magnitude smaller than the “typical values” normally 
adopted (§ 3.1). This reduced the effects of the dust grains, making our test model effectively 
a three-component plasma. We also neglected the source terms in the continuity equations, 
since chemistry was also neglected by Chernoff. The source terms for momentum transfer due 
to collisions between ions and neutrals were taken from DRD83 and Ciolek & Mouschovias 
(1993). 

For initial conditions, we used an adiabatic inflow model with a Mach number M = 9.75 
incident on a uniform initial state. The values for the density, temperature, and magnetic 
held strength of the inhow at z = 0 were determined by the jump conditions for a perpen¬ 
dicular, ideal (i.e., frozen-in magnetic hux) MHD J shock (e.g.. Priest 1984). The Alfven 
speed in the neutrals for the pre-shocked initial state is 6.17Cno- 

As the how expands into the computational region for t > 0, the initial J shock rapidly 
evolves into a C shock. The transient, short-term evolution during the transition from a 
J-type to C-type shock is similar to that which was reported by Toth (1994), and Smith & 
Mac Low (1997), including an early, rapid dihusion of the ions and magnetic held outward 
from the initially constricting J shock (with viz v^z in this phase of the evolution — for 
more discussion of transient phenomena in a representative model, see § 3.3.5 below). 

A steady C shock is soon established, with a velocity Vg = 12.6C'no. The Alfven Mach 
number M^a (= t's/'i^A.no) foi’ this shock is 2.04; this result is consistent with the analysis of 
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Chernoff (1987), who found that under these conditions (‘energy conserving’ flow), C shocks 
should exist for M^a < 2.8. Figure 8a shows the prohle of the neutral velocity (circles) 
and the ions (hlled diamonds) at a time after the steady state is attained. The smooth, 
supersonic-to-supersonic transition in the flow, characteristic of a C shock, is apparent in 
this hgure, and so is the magnetic precursor in the ion fluid that precedes the neutral flow 
at the head of the shock. 

The earlier studies of MHD shocks in partially ionized fluids were usually carried out in 
the frame of reference moving with the shock. Figure 8b displays the shock-frame velocities 
of the neutrals (un^,shk = Vnz — Vg] circles), and ions (uiz,shk = Viz — Vg] hlled diamonds), for 
the same time as in Figure 8a. As in previous examples, the velocities are normalized to 

C*nO- 

We may also quantitatively compare our numerical results to the detailed predictions of 
the phase-plane analysis of shocks in three-component plasmas by Chernoff (1987). Figure 
8 c displays the values (circles) of the phase variable r = |un 2 ,shk|/'ys as a function of g = 
I'i'i 2 ,shk|/'i's, for the numerical results presented in Figure 8b. The dashed lines in Figure 8c 
are the analytical values for r(g) calculated by using equation (35) of Chernoff (1987), for a 
shock with an Alfven Mach number M^a = 2.04. The essential equality of these two curves 
in this hgure indicates that the results of our multihuid numerical code are in accord with 
the analytical solution. 


3.3.5. Illustrative Solution with Grain Dynamics 

We now present an exemplary model that demonstrates some of the unique features 
of our time-dependent, multihuid numerical MHD code. As described above in § 2, this 
includes the ability to follow non-steady, transient phenomena, as well being able to follow 
simultaneously the dynamics of several diherent huid species. 

The illustrative model we present is an inhow J shock. The how for < 0, f < 0, is 
assumed to have a frozen-in magnetic held (so that the plasma, magnetic held, and neutrals 
move together), with the velocity and other quantities ai z = 0 determined by the jump 
conditions for a perpendicular shock. The resulting how inside the computational domain 
[0, Z] is then determined by solving the full set of governing equations, including the conti¬ 
nuity equations for all 9 huids (with source terms), and also the analytical cooling function 
An (§ 3.3). 

For specihcity, the pre-shock neutral gas has a uniform number density Uno = Pno/^n = 
1.03 X 10® cm“®, and temperature Tno = 10 K. The magnetic held of the uniform initial 
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state in this model is Bq = 600 /xG. The initial abundances of the various particle species 
are = 2.04 x 10"^°, = 3.26 x lO'^^ a;g-,o = 4.13 x XgO^o = 2.34 x lO'^^ 

Xg+^o = 3.42 X 10“^^, a;sg-_o = 9-94 x 10“^°, ^ggO g = 9.82 x 10“®, and Xgg+^o = 7.93 x 10“^'^. 
The density and magnetic held strength of this state are such that the Hall parameter for 
the small charged grains is Tgg = flsgT^g ^ = 25.4 3> 1. This means that the magnetic held 
in this model is ehectively frozen into the plasma of ions and also the small charged grains. 
The Hall parameter of the large grains is Tg = HgTg^ = 2.54 x 10“^ 1, so they are not 

attached to the magnetic held; this means that their dynamics are primarily governed by 
collisions with the neutrals. Finally, the inhowing J shock has M = 40, and M^a = 115, 
corresponding to a speed Vg = 9.8 km s“^. 

We hrst examine the model how on a length scale ~ 10^^ cm. On this length scale, 
and for these parameters, the shock travels as a typical J shock. Figures 9a-9c shows the 
prohle of various physical quantities at t = 5.5 x 10®s. Over these length and time scales, 
the coupling between the neutrals and the other species is rapid enough, compared to the 
characteristic how time of the shock, that the overall motion of the huids is essentially as a 
single, combined huid. 

In contrast, we now consider the same model on the much smaller scale 10^^ cm. 
This scale is now small enough that the ehect of drift of the magnetic held and the plasma 
with respect to the neutrals is discernible. So too, is the drift or lag of the large dust grains 
with respect to the neutrals because of inertial ehects. This is demonstrated in Figures 
10 a and 10c, which show the velocities of the neutrals (circles), ions and small charged 
grains (hlled diamonds), and the large dust grains (crosses) at times t = 1.32 x lO^s and 
t = 3.68 X lO^s, respectively. As can be seen from these hgures, these times are so early on 
that very little of the inhowing neutral gas has penetrated into the physical/computational 
region. However, there has been signihcant readjustment of the ions and the small charged 
grains, due to plasma drift. The dihusion of the charged particles is now resolved; the plasma 
and magnetic held have become ‘unlatched’ from the inhowing neutrals, as evident by the 
higher-velocity precursor in these hgures. 

Figure 106 shows the prohles of the neutral density (circles) and the magnetic held (hlled 
diamonds), normalized to the uniform reference values, at t = 1.32 x lO^s; Figure lOd shows 
the same thing at f = 3.69 x lO^s. Again, there is little penetration of the neutrals into 
the computational region during this early phase of the evolution. Because the magnetic 
gradient is much steeper at earlier times (being greatest at time t = 0, when the gradient 
is initially the discontinuity of the J shock), the velocity of the ions and the small charged 
grains is much larger at earlier times. This accelerates the plasma ahead of the shock, and 
gives rise to the precursor. The driving gradient decreases with time as the magnetic held 
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spreads out into the lower-density, lower-magnetic field strength undisturbed region ahead 
of the shock, and the velocity of the plasma correspondingly decreases from its maximum. 
From Figures 10a —d the effective signal or disturbance speed of the precursor at these early 
times is found to be ~ 10^ km s“^. This value is intermediate between the ion magnetosound 
speed and the magnetosound speed for the small charged grains, which are 1800 km s“^ and 
36 km s“^, respectively. 

Notably lagging behind the neutrals in Figures 10a and 10c are the large dust grains; 
only a few non-zero velocity data points for the large grains (crosses) appear in the lower 
left-hand corners of these figures. This is because of the inertia of the large grains, which 
is noticable on these scales. The distance over which the large grains are accelerated due to 
collisions with the neutrals is Ldrag(g) = — — Psdg/^nnn, where we have used the 

previously defined collision time scale, the relation nig = {in/3)pgttg^, and the fact that the 
relative drift speed in the shock is much greater than the sound speed. The neutral density 
in the shock front is ~ duno, and using = 2.33 amu for the assumed H 2 gas with a 10% 
abundance of He by number, we therefore have Ldrag(g) ~ 0.3 x lO^^cm. Hence, the lag of 
the large grains behind the neutrals on these length scales is understood for this model. 

The state of this model at t = 2.63 x lO^s is shown in Figures 11a and 116, and also 
at t = 5.5 X lO^s, in Figures 11c and lid. By these times the adjustment by plasma drift 
is approaching an asymptotic state, with an ion and small-charged grain mediated magnetic 
precursor still ahead of the shock. A particularly interesting feature in the magnetic field 
profile can be seen in Figures 116 and lid, which show that, because of the diffusion of the 
plasma with respect to the neutrals, magnetic field has “leaked forward” ahead of the shock. 
This decreases the magnetic field strength immediately behind the shock front relative 
to the far downstream (frozen-in) value ~ iB^. 

The lag of the large dust grains (crosses) behind the neutrals (circles) is again apparent 
in Figures 116 and lid. The distance over which the motion of the large grains is delayed 
with respect to the neutrals, because of large grain inertia, is consistent with, and explained 
by, the value for Ldrag(g)- 

Finally, we point out that the results presented in Figures 9a-9c are entirely consistent 
with the model results shown in Figures lOa-lOd and lla-lld. This is so because the small- 
scale structure has a maximum extent that is basically the size of the magnetic precursor, 
which has a width < 10^^ cm. However, examination of Figures 10 and 11 reveal that on 
scales greater than this, the values of the various physical quantities are the same as the 
upstream and downstream values for the J shock. Hence, on the large scale the smaller 
substructure is unresolved, and the model appears as a discontinuous MHD shock, such as 
that presented in Figure 9. 
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This model demonstrates our ability to follow the formation of multifluid MHD shocks, 
including short-lived phenomena. Detailed studies, for a wide variety of parameters and 
physical conditions, will be presented in later work. 


4. Summary 

In this paper — the hrst of a series — we have discussed the formulation of a time- 
dependent, numerical MHD code designed to model the formation and evolution of perpen¬ 
dicular MHD shocks in dense molecular gas. Chief among the various features of our code 
is that it is a multifluid code, incorporating the dynamics of neutral, ion, and six different 
types of grain fluids. The grain fluids consist of populations with two distinct radii, each 
with three different charge states. In addition to including for the effects of elastic collisions 
with gas particles, we also account for the transfer of mass and momentum between the grain 
species due to capture of charged particles on grain surfaces. We have presented the equa¬ 
tions required to model this physical system, and also described simplifying approximations 
adopted within the numerical code. 

To establish the accuracy of our numerical code, we have presented several test models 
and compared their results to exact, analytical solutions. In general, we found good agree¬ 
ment between the numerical and analytical solutions. This included a test of the formation 
of a C shock in a dust-free, weakly ionized gas, whose analytical solution was determined by 
Chernoff (1987). 

We also presented an illustrative model of a multifluid flow in which inertia of the large 
grains was shown to be important over a hnite scale of the flow. The possible effects that 
differential motion of fluids in shocks may have on the state of the neutral gas and other 
species will be investigated in more detail elsewhere. The effect of grains (including their 
inertia) on the speed and propagation of signals in dense gas will also be a topic of discussion 
in future work. Using this numerical code as a tool to aid in further research, we intend in 
subsequent papers to study the relevant physics and chemistry of the formation and evolution 
of various types of nonlinear MHD waves and shocks in molecular clouds, protostellar cores, 
and disks. 
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Table 1. 


Reaction 


Rate Coefficient (cm^ s 


cosmic ray + H 2 —^ m"*" + e a 


m+ + e —products 

m’*' + a —a"*" + m 

m+ + g~ —> m + g° 

m+ -|- gO —^ m + g+ 

m+ + sg“ —> m + sg° 
m"*" + sg° —> m + sg+ 


a+ 

+ 

e — 

a 

+ hu 

a+ 

+ 

g" - 

—a 

+ g° 

a+ 

+ 

gO _ 

a 

+ g+ 

a+ 

+ 

sg- - 

—a 

+ sg° 

a+ 

+ 

SgO _ 

—a 

+ sg+ 


e + gO —, g 

a+ + g“ —> a + g° 
m+ + g“ —> m + g° 
e + g+ —^ g° 

e + sg° —sg“ 

e + sg"*" —sg° 

g- + g+ —^ 2g° 

sg“ + sg+ —2sg° 


1.1 X 10-^(300/T) 

2.8 X 10"® 

b 

b 

b 

b 

2.45 X 10-12 (300/T)0-86 

b 

b 

b 

b 

b 

b 

b 

b 

b 

b 

b 

b 


^Ionization rate ( with C ~ 10 s ^ an adjustable parameter. 

I’Rate coefficient computed by integrating cross sections from Draine & Sutin (1987) over 
velocity distribution including relative drift. 
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Fig. la.— Density-dependent fractional abnndance of electrons (e), atomic ions (a"*"), molec- 
nlar ions (m"*"), small grains (sg) with radius Ogg = dA, and large grains (g) with radius 
Og = 0.1 /im, as predicted by the simplihed chemical network in Appendix A. The total abun¬ 
dances of large- and small grains in all charge states are Xg = 2 x 10“^^ and Xgg = 2 x 10“’^, 
respectively. 
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Fig. 2a.— Characteristic length scales for a steady shock with speed Vg = 10 km s plotted 
as a function of preshock density (see text). Small grains are present with Xgg = 2 x 10“^. 
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Fig. 3a.— Phase velocity of a small-amplitude wave propagating perpendicular to the mag¬ 
netic held in a cold ion-|-neutral plasma. Solid curves: exact dispersion relation. Open 
circles: effect of neglecting ion inertia. 












Fig. 4.— Nonmagnetic shock tube problem at time ti > 0. Initially the fluid is stationary, 
with Pn/Pno = 12, Tn/Tno = 1.25 for z/Z < 0.5, and Pn/Pno = 1, Fn/Tno = 1 for z/Z > 0.5. 
Circles: output from the numerical code described in this paper. Solid curves: exact solution, 
(a) Density, as a function of position z/Z. (b) Gas velocity, in units of the isothermal speed of 
sound Cno [= of the undisturbed gas. (c) Pressure, (d) Temperature. Since 

magnetic effects are not present in this particular model, all of the fluids move together as 
a single fluid because of effective collisional coupling. Note: for visual clarity, the curves 
showing the numerical results use only one out of every eight data points. 
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Fig. 5.— Same as in Figure 4, except at a later time t 2 > ti. By this time the forward¬ 
traveling shock wave and contact discontinuity have passed beyond z/Z = 1, and the 
backward-traveling rarefaction wave has extended to z/Z < 0. The absence of spurious 
numerical oscillations at z/Z = 0 and 1 is due to the non-reflecting condition implemented 
at the boundaries. 
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z/Z z/Z 


Fig. 6.— High velocity, adiabatic, nonmagnetic inflow test. At time t = 0 a shock with 
velocity Wnz = SOCno, Pn = 3.99pno, and = 835Tno enters the computational domain at 
z = 0. Shown is the model numerical solution (circles) at t > 0, by which time the shock has 
propagated into the computational domain, (a) Density, as a function of z/Z. (b) Velocity, 
(c) Pressure, (d) Temperature. The solid curves are the exact solution, which are shown for 


comparison. 
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Fig. 7.— Counterpropagating MHD rarefaction waves in a flux-frozen fluid. Initially, 
Pvl !PnO ^n/^nO Bx/B q 1, Unz/d^nO f fo^ z/Z < 1/^) aud V^iz/C^iQ 1 for 
z/Z > 1/2. Shown is the flow at a time t > 0, when the oppositely traveling rarefac¬ 
tion waves have evacuated a signihcant amount of mass and magnetic flux from the central 
portion of the physical region. Circles represent the numerical solution, which are compared 
with the MHD Riemann solution (solid curves) of Ryu & Jones (1995). Shown here are, nor¬ 
malized to the initial state variables, the prohles of the (a) density, {h) velocity, (c) pressure, 
and (d) magnetic held at this time. 
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Fig. 8.— An adiabatic MHD C shock in a weakly ionized gas, created by an inflow J shock 
with M = 9.75 incident on the boundary z = 0 at time t = 0. Model is shown after a steady 
C shock is attained, with speed Vs = 12.QC-no and an Alfven Mach number M^a = 2.04. 
(a) Velocities of the neutrals (circles) and the ions (filled diamonds). Characteristic of a C 
shock, there is a smooth, continuous transition from supersonic to subsonic flow speeds in the 
neutral fluid. A magnetic precursor in the ions runs ahead of the neutrals, {h) Velocities of 
the neutrals nnz,shk (circles) and ions Ui^^shk (Ailed diamonds) in the frame comoving with the 
shock. Shown are the negatives of the relative velocities nnz,shk = ivz — Vs) and Ui^^shk = {viz — 
Vs)i at the same time t as in (a), (c) Phase-plane variables r = |nn^,shk|/'ys and q = |niz,shk|/ns. 
Circles: numerical model results presented in (b). Dashed line: the theoretical prediction 
for a hydromagnetic C shock with = 2.04, calculated from equation (35) of Chernoff 
(1987). 
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Fig. 9.— Illustrative model with M = 40, MnA = 115. The model is shown at a time 
t = 5.5 X 10® s. (a) Velocity of the neutrals (circles), and ions and small charged grains (hlled- 
triangles). On these scales the velocities of the large grains are also the same as the neutrals, 
as are the small neutral grains, {h) Density of the neutrals, normalized to the upstream value 
Uno = 1-1 X 10® cm“®. (c) Magnetic held strength, normalized to Bq = 600 /xG. On these 
length scales, all of the huids propagate as a J shock. 
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Fig. 10.— Illustrative model, on much smaller length scales. The parameters for this model 
are the same as in the preceding hgure. (a) Velocities of the neutrals (circles), ions and small 
charged grains (hlled diamonds), and large dust grains (crosses) at time t = 1.32 x 10^ s. 
{h) Prohle of the density (circles, normalized to Uno) and magnetic held (hlled diamonds, in 
units of Bq) at the same time as in (a), (c) Same as in (a), except ed. t = 3.69 x 10^ s. (d) 
Same as {b), but at the same time as (c). At these early times, the magnetic held and the 
attached plasma species (ions, electrons, and small charged grains) rapidly dihuse ahead of 
the imposed hux-frozen how for z < 0 and t < 0, at a rate much faster than the inhowing 
nentrals and entrained large grains. 
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Fig. 11.— Ilustrative model shock, same as Figure 10 (including symbols), except at later 
times, (a) Velocities at t = 2.63 x lO^s. (b) Density and magnetic field structure, at same 
time as in (a), (c) Velocities at t = 5.50 x lO^s. (d) Density and magnetic field strength, 
same time as in (c). 







